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ABSTRACT 

In  the  conventional  atmospheric  turbulence  numerical  simulation  techniques  based  on  the  “split” -operator  method, 
the  turbulence -induced  refractive  index  inhomogeneities  are  represented  by  a  set  of  infinitely  narrow  (2D)  phase 
distorting  layers  (phase  screens).  These  conventional  phase  screens  cannot  represent  large-scale  refractive  index 
inhomogeneities  due  to  limitations  imposed  by  computational  grid  size.  For  this  reason,  this  commonly  used  model 
cannot  be  applied  for  computer  analysis  of  atmospheric  optical  systems  that  are  affected  by  the  presence  of  large- 
scale  turbulence  eddies.  Among  these  systems  are  coherent  imaging  ladars,  optical  vibrometers  and  interferometers. 
In  the  classical  Kolmogorov  turbulence  theory,  the  impact  of  the  large-scale  turbulence  eddies  are  associated  with 
the  turbulence  outer  scale.  Contrary  to  the  conventional  approach,  in  the  computer  simulation  technique  introduced 
here  the  turbulence -induced  refractive  index  inhomogeneities  are  represented  by  a  set  of  large-scale  phase  distorting 
screens  that  account  for  refractive  index  inhomogeneities  which  extend  beyond  the  numerical  grid  correlation 
length.  The  results  are  applied  to  the  analysis  of  piston  phase  fluctuations  for  the  cases  when  the  turbulence  outer 
scale  significantly  exceeds  the  receiver  aperture  size.  We  also  analyze  the  piston  phase  fluctuations  in  deep 
turbulence  conditions  in  presence  of  phase  singularities  (phase  cuts  and  branch  points),  and  show  that  the 
conventional  definition  of  piston  phase  cannot  be  applied  for  this  case.  We  introduce  a  more  general  definition  of 
piston  phase  which  is  useful  for  analysis  in  deep  turbulence  conditions. 

1.  INTRODUCTION 

The  efficiency  of  various  optical  systems  based  on  coherent  detection  techniques  heavily  depends  on  the  accurate 
prediction  of  statistical  and  temporal  characteristics  of  the  aperture  averaged  (piston)  phase  [1].  The  rigorous 
estimation  of  the  piston  phase  represents  a  quite  challenging  task,  since  the  piston  phase  depends  strongly  on  large- 
scale  turbulent  eddies  with  sizes  on  the  order  of  or  exceeding  the  receive  aperture  diameter  D  [2].  Although  the 
smaller  size  turbulent  eddies  can  also  affect  piston  phase  values,  their  impact  is  less  profound  due  to  aperture 
averaging.  Thus,  accurate  estimation  of  piston  phase  cannot  be  performed  without  consideration  of  the  largest  size 
refractive  index  fluctuations.  Within  the  framework  of  the  fully-developed  Kolmogorov  turbulence  model,  the 
largest  scale  eddies  are  described  by  the  turbulence  outer  scale  L0,  whose  size,  ranging  from  a  few  to  hundreds  of 
meters,  can  significantly  exceed  the  receiver  aperture  diameter  D. 

Predictive  modeling  of  optical  wave  propagation  and  piston  phase  estimation  is  currently  performed  using  the 
conventional  representation  of  Kolmogorov’s  turbulence  by  a  set  of  statistically  independent ,  infinitely  thin  two- 
dimensional  phase  screens,  which  are  spatially  hounded  inside  a  numerical  grid  domain.  This  conventional  approach 
creates  significant  problems  in  the  analysis  of  optical  systems  (including  coherent  detection  systems)  whose 
performance  depends  on  large-scale  turbulent  eddies.  In  the  conventional  numerical  simulation  techniques  used  for 
phase  screen  generation,  the  turbulence  outer  scale  L0  is  associated  with  the  numerical  grid  size  and  exceeds  D  only 
by  a  few  fold  (typically  from  four  to  eight  fold  only).  Because  of  practical  considerations  related  to  acceptable 
duration  of  the  computation  time  and  available  computer  memory,  in  most  predictive  wave -optics  simulations  the 
grid  size  typically  doesn’t  exceed  1024x1024  (more  commonly  512x512)  and  accordingly,  the  corresponding 
analysis  is  limited  by  atmospheric  conditions  for  which  the  outer  scale  L0  is  on  the  order  of  the  receiver  aperture 
diameter  D.  Clearly,  piston  phase  estimations  based  on  wave  optics  simulations  that  utilize  such  conventional  phase 
screen  generation  techniques,  may  not  fully  account  for  the  impact  of  large-scale  turbulent  eddies  and  this  may  lead 
to  significant  errors. 

Another  serious  problem  arises  in  modeling  of  the  temporal  dynamics  of  the  piston  phase  in  tracking  of  moving 
objects.  In  this  application,  wave  propagation  occurs  under  conditions  of  continuously  changing  optical  axis 
direction  and  potential  crossing  of  several  atmospheric  layers,  which  may  have  different  turbulence  properties  over 
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the  propagation  path.  All  these  propagation  scenarios  require  consideration  of  a  large  region  of  the  atmosphere, 
which  may  need  significantly  large  almost  impractical  numerical  grid  size.  In  this  paper  we  address  these  problems 
by  applying  for  piston  phase  analysis  the  recently  developed  infinitively  long  phase  screen  (ILPS)  technique  [3].  In 
the  following  section  we  briefly  describe  the  ILPS  technique.  Analysis  of  piston  phase  fluctuation  variance  is 
presented  in  section  3. 

2.  GENERATION  OF  PHASE  SCREENS  WITH  ARBITRARY  LARGE  TURBULENCE  OUTER  SCALE 

In  the  conventional  phase  screen  generation  technique,  a  random  function  S  =  S ( m,n )  ,  defined  inside  a  square 
grid,  is  computed  using  Fourier  transform  of  the  product  of  refractive  index  fluctuation  power  spectrum  function 
O (&,/)  and  the  delta-correlated  on  the  Ax  A  grid  complex  random  function  rj(kj)  with  uniform  probability 
distribution  inside  [0,  1]  interval,  as: 

\ 

S(m,n)  =  (1) 

where  (m,n)  ,  m,n  e  {l,...,  A}  and  (&,/)  ,&, l  e  {l,...,  A}  )  are  integer  numbers  corresponding  to  numerical  grid 
pixels  and  FT  is  the  Fourier  transform  operator.  For  statistically  independent  real  and  imaginary  components  of 
function  r/(k,l)  the  complex  function  S  is  comprised  of  two  statistically  independent  (component)  phase  screens. 
These  functions  are  periodic  with  the  numerical  grid  size  period.  The  examples  of  the  conventional  phase  screens 
are  shown  in  Fig.  1. 


Fig.  1.  Grey  scale  representation  of  the  generated  on  512x512  grid  random  phase  screen  realizations  with  different  power 

spectrums:  Kolmogorov  (left)  and  Tatarskii  (right). 

In  the  next  step,  consider  computer  generation  of  a  phase  screen  S  that  is  defined  inside  an  extended  in  one 
direction  grid  of  size  JN  x  A ,  where  J  is  the  grid  extension  factor  (J  =  2, 3, ... ).  As  shown  in  [4]  such  phase  screen 

can  be  obtained  by  summarizing  J  conventional  random  phase  screen  functions  Sj{m,n)  ( j  e  {l, . ..,/})  as 
defined  by  expression  (1)  with  exponential  weighting  factors: 


j  2njm 

S(m,n)  =  Ye  JN  S.{m,n),  (2) 

7=1 

where  m  e  {l,..,/A}  and  n  e  {l,..,  A}  ,  and  i2  =  -1 .  Note  that  the  obtained  new  phase  screen  S  is  a  periodic 
function  with  the  period  JN  along  x  and  A  along  y  directions.  In  the  coordinate  space,  the  extended  grid 

corresponds  to  LxLy  area,  where  Lx  -  dxJN  and  Ly  -  dyN  and  dx  and  dy  are  pixel  lengths  in  physics  units. 
Examples  of  phase  screens  obtained  on  an  extended  grid  (extended  phase  screens)  are  shown  in  Fig.  2. 


Fig.  2.  Gray  scale  representation  of  the  generated  extended  phase  screens  corresponding  to  Kolmogorov  (top)  and  Tatarskii 
(bottom)  power  spectra  on  the  numerical  grid  containing  4*512x512  pixels  (the  grid  extension  factor  J= 4). 
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Generation  of  phase  screens  with  arbitrary  large  turbulence  outer  scale  L0  can  now  be  performed  by  utilizing  the 
extended  phase  screen  technique  described  above.  The  idea  is  based  on  partitioning  of  the  spectral  domain  of  the 
refractive  index  power  spectrum  function  ®  with  a  predefined  outer  scale  value  L0  onto  H  >  1  ( H  is  an  integer 

number)  spectral  sub-regions  and  generation  of  the  extended  phase  screens  Sh  (r)  ,  where  h  e  {l,...,//}  ,  with  the 

power  spectrum  functions  ®/7  (k)  that  coincide  with  ®  inside  these  sub-regions  and  equal  to  zero  otherwise.  Here 
k  is  a  vector  in  the  spectral  domain.  Consider  as  an  example  partitioning  of  the  spectral  domain  onto  sub-regions  in 
the  form  of  concentric  annuli  centered  at  the  coordinate  origin.  As  shown  in  [3],  phase  screens 

corresponding  to  the  power  spectrum  ®  with  outer  scale  L0  can  be  obtained  using  summation  of  the  extended  phase 
screens  Sh  (r)  corresponding  to  the  power  spectrum  functions  ®/7  (k)  : 

S(r)  =  l>*(r)  (3) 

h= 1 

The  extended  phase  screens  Sh  (r)  are  defined  inside  coordinate  domain  Lx  -  dxJN  >  L0  ,  L  -  dyN  and  are 
generated  using  the  described  above  technique.  An  example  of  phase  screen  with  Tatarskii  power  spectrum  and 
L0  =  1L  is  presented  in  Fig.  3.  Note  that  fusion  of  the  extended  phase  screens  Sh  (r)  that  are  generated  inside 

coordinate  domains  with  different  pixel  sizes  dxh  and  dyh ,  requires  additional  computation  (approximation  and 
filtering)  to  unify  grid  pixel  sizes  prior  to  the  fusion  [3]. 


i 

Fig.  3.  Gray  scale  representation  of  phase  screen  corresponding  to  Tatarskii  power  spectrum  with  the  outer  scale  L0  that  is  seven 
fold  longer  than  the  phase  screen  width.  The  phase  screen  is  generated  using  numerical  grid  with  8*512x512  (J=8)  pixels. 

In  any  numerical  simulations  of  atmospheric  optical  systems,  temporal  dynamics  of  the  atmospheric  turbulence  may 
play  an  important  role.  Conventional  approach  to  include  turbulence  induced  temporal  variations  in  system 
performance  is  based  on  Taylor’s  hypothesis  of  “frozen”  refractive  index  inhomogeneities  that  are  moved  as  a  whole 
with  a  wind  velocity  v  .  In  the  split  operator  technique  described,  the  wind  induced  temporal  variations  inside  optical 
system  aperture  are  simulated  by  introducing  a  set  lateral  shifts  of  phase  screens  over  distance  dx  ,  that  corresponds 

to  the  grid  pixel  size  in  the  direction  orthogonal  to  the  system  optical.  The  characteristic  time  TG  for  the  complete 
update  of  the  extended  phase  screen  realization  that  is  defined  inside  coordinate  domain  Lx  -  dxJN  >  L0,Ly=dyN 
equals  to  TG  =  Lx/vx  <  L0/vx,  where  is  x-  projection  of  wind  velocity  vector  v.  For  simplicity  we  assume 
that  vy  =  0  .  Thus  the  conventional  approach  of  phase  screens  generation  allows  only  analysis  of  temporal  dynamics 

over  time  T  <TG  that  is  associated  with  impact  of  a  single  turbulence  eddy  of  size  L0.  This  restriction  significantly 
limits  ability  to  compute  time-averaged  characteristics  that  are  dependent  on  the  turbulence  outer  scale,  such  as 
piston  phase,  and  hence  require  computer  simulation  of  temporal  dynamics  over  the  time  T  <TG  that  significantly 

exceeds  the  time  of  the  extended  phase  screen  update  time.  As  shown  in  [3],  this  problem  can  be  resolved  by 
applying  the  infinitively-long  (IL)  phase  screen  technique. 

Assume  a  set  of  the  extended  phase  screens  =  S"'  ( m,n ) ,  where  j  e  {l,...,//L}  .  Here  JIL  =  pJ  and  p  >2  is  an 

integer.  Each  of  these  extended  phase  screens  are  generated  using  the  described  above  technique  and  are  defined 
inside  the  extended  grid  with  JN  x  N  pixels.  Consider  the  following  defining  the  infinitely  long  phase  screen 
expression 
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(4) 


JIL 

SIL  ( m,n )  =  (m,n)S^  (m  -  jM,n )  , 

7=0 

where  y/  ( m,n )  =  cos (nm  /  / /ZA)  ,  M  =  J ILN  /  2 ,  m  >  1  and  n  e  {l,...,  TV} .  Note  that  the  S/L  =  *S/L  ( m,n )  is 
defined  inside  an  infinite  in  x  -direction  grid.  As  shown  in  [4],  the  random  function  SIL  defined  by  (4)  has 
approximately  identical  statistical  properties  as  random  functions  .  Note  that  in  each  grid  point  ( m,n ) 

function  SIL  ( m,n )  is  composed  of  only  two  statistically  independent  extended  phase  screen  functions  multiplied  by 
the  weighting  factors  y/  ( m,n )  and  y/j+l  ( m,n )  :  random  function  realization  S ^  and  the  shifted  over  M  pixels  in 
v -direction  random  function  S^+1  as  illustrated  in  Fig.  4.  This  means  that  for  computer  simulation  of  temporal 
dynamics  over  an  arbitrary  long  time  T  >TG  using  sequential  shifts  of  phase  screens,  there  is  no  need  for 

generation  and  keeping  in  the  computer  memory  of  large  number  J IL  of  random  functions  S".  ,  ( j  e  {l, ...,  J IL  } ).  In 
fact,  the  same  goal  can  be  achieved  by  generation  and  saving  in  computer  memory  only  two  extended  phase  screens 
( and  S^+l ).  The  weighted  combination  of  these  functions  5  a  as  described  by  (4)  can  be  used  until  the  lateral 

shift  of  the  phase  screen  approaches  the  extended  grid  boundary.  At  this  point  an  additional  the  extended  phase 
screen  S^+2  should  be  generated.  The  temporal  dynamics  inside  the  next  time  interval  of  duration  TQ  is  now  defines 

by  phase  screen  SIL  that  contains  weighted  combination  of  functions  S"+l  andS^2 .  This  process  of  regeneration  of 

extended  phase  screens  and  their  fusion  can  be  repeated  as  many  times  as  required  for  analyzing  an  optical  system 
dynamics  over  arbitrary  long  time. 


Extended  phase  screen  SIL 


Extended  phase  screen  S"'+i 


Section  of  phase  screen  SIL 


Fig.  4.  Graphical  illustration  of  an  infinitely  long  phase  screen  generation  technique. 

3.  STATISTICAL  CHARACTERISTICS  OF  PISTON  PHASE:  NUMERICAL  SIMULATION  RESULTS 


3.1  Pupil-plane  phase  screen  model.  In  this  section  we  consider  statistical  characteristics  of  piston  phase  measured 
at  a  receiver  telescope  of  diameter  D  which  are  obtained  through  numerical  simulation  of  plane  wave  propagation 

through  a  moving  with  wind  velocity  v  =  (vx,0)  atmospheric  turbulence.  The  commonly  used  definition  of  piston 
phase  is  given  by  the  following  expression: 
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(5) 


(p  {t)  =  —  [  (p(r,t)d2r 
Sa 

where  <p(r,t)is  the  instantaneous  phase  function  defined  inside  the  receiver  telescope  aperture  area  Sa.  First 

assume  that  the  impact  of  atmospheric  turbulence  can  be  described  by  a  single  thin  phase  screen  that  is  located  at  the 
telescope  pupil  (pupil-plane  phase  screens).  This  model  is  commonly  used  for  analysis  of  optical  propagation  in  the 
Earth  atmosphere  at  high  elevation  angles.  In  the  numerical  simulations  of  piston  phase  dynamics,  we  used  pupil- 
plane  phase  screens  with  Tatarskii  refractive  index  fluctuation  power  spectrum  with  a  fixed  parameter  Cn 2  and 
different  values  of  the  turbulence  outer  scales  L0  ranging  from  L0  =  D  to  L0  =10 D.  The  numerical  simulations  were 
performed  for  receiver  aperture  of  diameters  D= 3.6  m  and  D=  1.0  m.  We  also  assumed  that  all  temporal  changes  in 
the  piston  phase  are  due  to  lateral  translation  of  the  phase  screen  with  a  constant  wind  speed  v=(vx=1.0  m/sec,  0).  In 
the  numerical  generation  of  phase  screens  the  grid  extension  factor  J  in  (2)  was  varied  depending  on  the  chosen 
value  of  the  outer  scale  value  L0.  We  simulated  piston  phase  dynamics  in  T  =24  min  time  intervals.  This  time 
interval  corresponds  to  translation  through  the  receiver  aperture  D  of  the  pupil  plane  turbulent  screen  of  length 

400 D.  The  total  length  of  effective  phase  screen  was  controlled  by  setting  parameter  J IL  in  (4). 

Examples  of  piston  phase  dynamics  for  two  different  values  of  the  outer  scale  are  presented  in  Fig.  5.  The  results 
suggest  that  increase  of  L0  leads  to  rapid  increase  of  piston  phase  fluctuations.  Using  the  piston  phase  temporal 
evolution  curves  as  in  Fig.  5,  we  estimated  the  standard  deviation  of  piston  phase  fluctuations 

<r  =  (^(t)  —  (^(t))J  \  .  In  computation  of  a ,  we  substituted  ensemble  averaging  by  time  averaging  over  24 

min  observation  time.  The  results  of  a  estimation  are  shown  in  Fig.  6  for  two  different  receiver  aperture  diameters. 
The  standard  deviation  of  piston  phase  fluctuations  rapidly  increases  with  the  outer  scale  increase  reaching  2tt  value 
at  L0  «  36m  in  the  example  considered.  At  the  same  time  cr  only  weakly  depends  on  the  receiver  aperture  size 
(compare  two  curves  corresponding  to  D= 3.6  m  and  D=  1.0  m  in  Fig.  6). 

m  * 
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—71  - 

-271  - 

-3  71  I - 1 - » - 1 - 1 - 1 - ( - 1 - 

0  3  6  9  12  15  18  21  24 
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Fig.  5.  Temporal  evolution  of  piston  phase  for  L0  =D  (left)  and  L0  =5D  (right)  for  Cn 2  =  10'15  m'2/3,  and  D= 3.6  m. 


Lq/Dq 

Fig.  6.  Piston  phase  standard  deviation  a  vs  turbulence  outer  scale  L0  for  Cn 2  =  10~15  m'2/3  . 

Consider  now  the  dynamics  of  piston  phase  for  slant  propagation  paths.  In  these  propagation  scenarios,  phase 
function  (p{r,t )  of  entering  receiver  aperture  optical  field  is  characterized  by  the  presence  of  intensity  scintillations, 
271-phase  cuts  and  branch  points.  The  first  question  to  ask  is  how  to  define  piston  phase  for  slant  propagation 
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conditions  by  taking  into  account  sensitivity  of  piston  phase  value  to  geometry  of  27i-phase  cuts  and  branch  points. 
Indeed  in  presence  of  these  phase  singularities,  the  commonly  used  definition  of  piston  phase  (5)  doesn’t  work  since 

phase  function  q>(r,t )  is  not  uniquely  defined.  Geometry  of  27i-phase  cuts  in  phase  function  (p{r,t )  can  be 

changed  without  having  any  impact  on  optical  field  complex  amplitude  that  is  proportional  to  exp^i0>(r,f)J  but  at 

the  same  time  this  change  affects  the  piston  phase  since  the  aperture  averaged  value  of  phase  function  in  (5)  depends 
on  27i-phase  cuts  geometry.  To  address  this  problem,  we  introduce  here  a  different  definition  of  piston  phase  that  is 
insensitive  to  presence  of  phase  singularities.  We  define  piston  phase  7px  (t)  as  solution  of  the  following  equation: 

|  sin(^(r,/)  -  cpl  (f))d2r  =  0 ,  (6) 

In  this  definition,  one  can  add  In  value  to  at  any  point  r  inside  receiver  aperture  without  changing  piston 
phase  cp]  (/)  .  The  solution  of  (6)  can  be  represented  in  the  form: 

sin(^(r,^))  dr 

7 - T  ’  (7) 

|  cos(<£?(r,£))  d2r 

Sa 

Note  that  although  the  phase  defined  by  (6)  is  insensitive  to  2n  phase  cuts  (jumps),  its  value  still  is  not  uniquely 

defined  because  function  tan-1  in  (7)  is  a  multivalued  function  with  infinite  number  of  separated  by  n  branches. 
The  presence  of  branch  points  can  result  in  uncertainty  in  selection  of  one  or  another  branch  of  this  function.  A 
special  technique  based  on  analysis  of  phase  function  and  location  of  phase  singularities  was  applied  to  minimize 
impact  of  2n  phase  cuts  and  branch  points  on  the  piston  phase  computation.  The  method  of  piston  phase 
computation  based  on  (7)  was  compared  with  the  corresponding  computations  using  conventional  piston  phase 

function  (p  (r,f)  .  The  comparative  analysis  was  performed  for  D=0.5  m  and  the  outer  scale  values  ranging  from 

L0=D  toL0=10£>.  As  illustrated  in  Fig.  7,  the  obtained  values  of  standard  deviations  for  the  piston  phases  cp  (f )  and 

(px  (t)  coincide  with  good  accuracy  which  supports  legitimacy  of  the  introduced  piston  phase  definition  (7). 

Consider  now  results  of  analysis  of  piston  phase  as  defined  by  (7)  in  deep  turbulence  conditions  characterized  by 
strong  scintillations.  In  the  numerical  simulations  we  used  three  statistically  independent  extended  phase  screens. 
The  first  phase  screen  was  located  in  the  receiver  telescope  pupil  plane  at  z=Zi=0  and  the  second  and  third  phase 
screens  at  distances  z2=  0.025  kD2  and  z3=  0.05  kD 2  from  the  receiver  plane.  All  three  phase  screens  were  literary 
translated  to  account  for  wind  speed  y=(vx=1.0  m/sec,  0). 
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Fig.  7.  Dependence  of  piston  phase  standard  deviation  a  on  turbulence  outer  scale  L0  for  the  case  of  pupil  plane  phase  screen 
(Tatarskii  power  spectrum,  Cn 2  =  10'15  m'2/3  and  D0=0.5  m).  The  top  curve  corresponds  to  the  conventional  piston  phase 
definition  (5)  and  the  bottom  line  to  the  piston  phase  definition  (7)  which  is  introduced  for  piston  phase  analysis  in  deep 
turbulence  conditions. 


<P,  (t)  =  tan 
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The  results  of  piston  phase  dynamics  over  24  min  time  interval  are  presented  in  Fig.  8  for  two  different  values  of  the 
turbulence  outer  scale. 


Fig.  8.  Temporal  evolution  of  piston  phase  in  deep  turbulence  conditions  that  are  generated  using  three  distant  phase  screens  for 
L0  =5D  (left)  and  L0=\OD  (right)  for  Cn2  =  10' 15  m'2/3,  and  D=0.5  m. 


Fig.  9.  Piston  phase  standard  deviation  a  vs  turbulence  outer  scale  L0  for  pupil  plane  and  distant  phase  screens 

(Cn2  =  10'15  m'2/3  and  D0= 0.5  m). 


Comparison  of  piston  phase  fluctuations  for  a  single  phase  screen  located  either  at  the  pupil  plane  or  shifted  a 
distance  z3=  0.05  kD 2  from  the  receiver  pupil  (distant  phase  screen)  are  presented  in  Fig.  9,  where  the  standard 
deviation  of  the  piston  phase  is  shown  as  a  function  of  the  outer  scale.  The  results  show  that  for  L0  <  5 D 
characteristic  range  of  piston  phase  fluctuations  as  measured  by  its  standard  deviation  is  nearly  identical  for  both 
pupil  and  remotely  located  phase  screens.  With  increase  of  the  outer  scale,  the  piston  phase  fluctuations  are  smaller 
for  the  remotely  located  phase  screen  (deep  turbulence). 
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